# Figure 2
rm(list=ls())
gc()

# Load RCPs
## First column is temperature time series from deterministic EBM, initialised at 1C in 2020. Next 10,000 columns are temperature time series from stochastic EBM. See Matlab script "figure2.m" for further details.
path <- "~/Documents/Projects/IAM_Stochastic/" # Directory where output from figure2.m is stored.
files <- list.files(path=path, pattern="output.txt")
d <- data.frame()
for (i in 1:length(files)){
	ds <- read.table(paste(path, files[i], sep=""), sep=",", header=FALSE)
	ds$rcp <- files[i]
	d <- rbind(d, ds)
}

# Associate RGB values for IPCC colour scheme with RCPs
d$r <- NA
d$r[grep("26",d$rcp)] <- 0
d$r[grep("45",d$rcp)] <- 121
d$r[grep("60",d$rcp)] <- 255
d$r[grep("85",d$rcp)] <- 255
d$g <- NA
d$g[grep("26",d$rcp)] <- 0
d$g[grep("45",d$rcp)] <- 188
d$g[grep("60",d$rcp)] <- 130
d$g[grep("85",d$rcp)] <- 0
d$b <- NA
d$b[grep("26",d$rcp)] <- 255
d$b[grep("45",d$rcp)] <- 255
d$b[grep("60",d$rcp)] <- 45
d$b[grep("85",d$rcp)] <- 0

# Ensemble size
nn = ncol(d) - 6
names(d) <- c("year","det",paste("X",seq(1:nn),sep=""), "rcp","r","g","b")
d$year <- d$year + 2019

# Length of time series
n=length(unique(d$year))

# Economic assumptions 
gdp = 80*10^12 # initial aggregate consumption (80 trillion US dollars)
g = 0.019 # annual growth rate of undisturbed aggregate consumption
eta = c(0,1.45) # marginal utility of consumption
rho = c(0.04255,0.015) # pure rate of time preference (1.5% in DICE-2016)

s <- 1 # Index for discounting assumptions. When s=1, r = 4.255 + 0 * 0.019 = 4.255.

# Demographic assumptions
pop = 7.5 # billion people
pop.asympt = 11.5 # Asymptotic population, in billion
param = 0.134/5 # DICE-2016 "Population 2050 parameter"
pop.ts <- rep(NA,n)
pop.ts[1] <- pop
for (i in 2:n) { pop.ts[i] <- pop.ts[i-1]*(pop.asympt/pop.ts[i-1])^param}
pop.ts <- matrix(rep(pop.ts,nn+1),ncol=nn+1,nrow=n)
pop.ts <- pop.ts*10^9 # Population

# Undisturbed output per capita
gdp.ts <- gdp*(1+g)^seq(0,n-1,1)
gdp.ts <- matrix(rep(gdp.ts,nn+1),ncol=nn+1,nrow=n)
gdp.ts <- (gdp.ts/pop.ts)/1000 # undisturbed consumption per capita expressed in thousands of dollars.

# Discount factors
df.2 <- 1/((1+rho[s])^seq(0,n-1,1))
df.2 <- matrix(rep(df.2,nn+1),ncol=nn+1,nrow=n)

x <- data.frame()
for (j in files) {
	t <- d[d$rcp==j,]
	t <- t[,!(names(t) %in% c("year","rcp","r","g","b"))]
	# Isoelastic utility
	x.Weitzman <- colSums((((gdp.ts*(1/(1+((t)/20.46)^2 + ((t)/6.081)^6)))^(1-eta[s]))/(1-eta[s])) * df.2 * pop.ts)
	x <- rbind(x, x.Weitzman)
}
x <- t(x)
colnames(x) <- files

## Scaling factors for plotting
UofC <- ((gdp.ts[1,1])^(1-eta[s]))/(1-eta[s]) # Utility of undisturbed per capita consumption today (in thousands of dollars)
scalar <- UofC-(((gdp.ts[1,1]-0.001)^(1-eta[s]))/(1-eta[s])) # Utility value of the marginal dollar
UofC <- sum((((gdp.ts[,1])^(1-eta[s]))/(1-eta[s]))*pop.ts[,1]*df.2[,1]) # NPV of utility of undisturbed aggregate consumption


# # Unused diagnostic plot: We simulated temperature ensembles of different sizes (using "figure2.m") to confirm that the estimates of economic damages converge. This code can be used to plot the estimated 5th and 95th percentiles of the damage distribution as a function of ensemble size.
# j = rev(files)[1]
# k = 2
# l = 1000
# y <- data.frame()
# while ((k+l) < nn) {
# 	dens <- density((UofC - x[k:(k+l),grep(j,colnames(x))])/scalar)
# 	y <- rbind(y, c(l, max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])))
# 	k <- k+l+1
# 	l <- l+1000
# }
# names(y) <- c("n","x5","x95")
# plot(y$n,y$x95,type="l", ylim=c(min(y$x5), max(y$x95)))
# lines(y$n,y$x5)

# Estimate risk premium
s <- 2 # Index for discounting assumptions. When s=2,  r = 1.5 + 1.45 * 1.9 = 4.255.

# Discount factors
df.1 <- 1/((1+rho[s])^seq(0,n-1,1))
df.1 <- matrix(rep(df.1,nn+1),ncol=nn+1,nrow=n)

x1 <- data.frame()
for (j in files) {
	t <- d[d$rcp==j,]
	t <- t[,!(names(t) %in% c("year","rcp","r","g","b"))]
	# Isoelastic utility
	x.Weitzman <- colSums((((gdp.ts*(1/(1+((t)/20.46)^2 + ((t)/6.081)^6)))^(1-eta[s]))/(1-eta[s])) * df.1 * pop.ts)
	x1 <- rbind(x1, x.Weitzman)
}
x1 <- t(x1)
colnames(x1) <- files

## Scaling factors
UofC1 <- ((gdp.ts[1,1])^(1-eta[s]))/(1-eta[s]) # Utility of undisturbed per capita consumption today (in thousands of dollars)
scalar1 <- UofC1-(((gdp.ts[1,1]-0.001)^(1-eta[s]))/(1-eta[s])) # Utility value of the marginal dollar
UofC1 <- sum((((gdp.ts[,1])^(1-eta[s]))/(1-eta[s]))*pop.ts[,1]*df.1[,1]) # NPV of utility of undisturbed aggregate consumption

j = rev(files)[1]
((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1) / 10^12 # Risk premium in present value dollars
((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1)/gdp # Risk premium as a share of current global output

j = rev(files)[2]
((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1) / 10^12 # Risk premium in present value dollars
((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1)/gdp # Risk premium as a share of current global output

j = rev(files)[3]
((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1) / 10^12 # Risk premium in present value dollars
((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1)/gdp # Risk premium as a share of current global output

j = rev(files)[4]
((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1) / 10^12 # Risk premium in present value dollars
((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x),grep(j,colnames(x1))]))/scalar1)/gdp # Risk premium as a share of current global output

# Plot
quartz(width = 10, height = 5, type = "pdf", "Figure4", file = "Figure4.pdf")
par (fig=c(0,1,0,1), # Figure region in the device display region (x1,x2,y1,y2)
	omi=c(.8,0,0,.5), # global margins in inches (bottom, left, top, right)
	mai=c(0,.8,0.4,.3)) # subplot margins in inches (bottom, left, top, right)
	layout(matrix(c(1,1,1,1,2,3,4,5), 4, 2, byrow = FALSE), widths=c(4,3), heights=c(3,3,3,3))

# Panel 1: Temperature trajectories
plot(d$year[d$rcp==files[1]],d$det[d$rcp==files[1]], type="n", xlim=c(2020,2200), ylim=c(0,10), bty="L", ylab="", xaxs="i",yaxs="i",yaxt="n",main="Global mean temperature anomaly")
mtext(side=1, "Year", outer=T, line=2, cex=.8, at=0.3)
mtext(side=3, expression(Delta * "T (˚C)" ), outer=T, line=-2.5, cex=.8, at=0.07)
axis(side=2,las=2)
for (j in files) {
	lines(d$year[d$rcp==j],d[d$rcp==j,i], col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 255))
	lines(d$year[d$rcp==j],d$det[d$rcp==j])
	mtext(paste(" ", substr(j,1,4), ".", substr(j,5,5), sep=""), side = 4, las=2, line = 0, outer = FALSE, at = d$det[d$rcp==j & d$year==2200],  adj = NA, padj = NA, cex = NA, col = NA, font = NA)
}

# Panels 2-5: Damage distributions
j = rev(files)[1]
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="Damages",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(min((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar),max((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar))) # xlog=T,xlim=c(10,6*10^14)
mtext(side=1, "Trillion USD", outer=T, line=2, cex=.8, at=0.81)

polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±5% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*max((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar), y=0.9*max(dens$y), labels="Risk Premium", cex=1.2, pos=2)
text(x=0.99*max((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar), y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*max((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar), y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)


j = rev(files)[2]
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(min((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar),max((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)))
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±10% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*max((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar), y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*max((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar), y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)


j = rev(files)[3]
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(min((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar),max((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)))
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*max((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar), y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*max((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar), y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)


j = rev(files)[4]
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(min((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar),max((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)))
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±10% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*max((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar), y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*max((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar), y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)

dev.off()

